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Abstract: A new class of Marginal Structural Models (MSMs), History- 
Restricted MSMs (HRMSMs), was recently introduced for longitudinal data 
for the purpose of defining causal parameters which may often be better 
suited for public health research or at least more practicable than MSMs 
B 0) • HRMSMs allow investigators to analyze the causal effect of a treat- 
ment on an outcome based on a fixed, shorter and user-specified history 
of exposure compared to MSMs. By default, the latter represent the treat- 
ment causal effect of interest based on a treatment history defined by the 
treatments assigned between the study's start and outcome collection. We 
lay out in this article the formal statistical framework behind HRMSMs. 
Beyond allowing a more flexible causal analysis, HRMSMs improve com- 
putational tractability and mitigate statistical power concerns when de- 
signing longitudinal studies. We also develop three consistent estimators 
of HRMSM parameters under sufficient model assumptions: the Inverse 
Probability of Treatment Weighted (IPTW), G-computation and Double 
Robust (DR) estimators. In addition, we show that the assumptions com- 
monly adopted for identification and consistent estimation of MSM param- 
eters (existence of counterfactuals, consistency, time-ordering and sequen- 
tial randomization assumptions) also lead to identification and consistent 
estimation of HRMSM parameters. 

Keywords and phrases: causal inference, counterfactual, marginal struc- 
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1. Motivation 

Longitudinal epidemiological studies are increasingly becoming more interested 
in the time-dependent effects of various exposures on human health outcomes. 
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This is particularly true for studies concerned with the effects of chronic ex- 
posure to ambient air pollutants. Cohort studies with multiple time-specific 
estimates of exposure have been emphasized by the U.S. Environmental Pro- 
tection Agency (EPA) as the preferred study design to address these issues (|l|). 
Examples of such concerns can be found for other health relevant exposures, 
such as time-specific patterns of physical activity on cardiovascular outcomes 
and obesity (j22l ). 

Typically, cohort studies collect data at regular time intervals for all members 
of the cohort. In practice, each collection time represents a "window" of time over 
which data are collected. Information on the exposure of interest, also referred 
to as "treatment of interest" , and other relevant covariates are obtained for the 
interval between successive data collection time points. 

Currently, the principal tools used by epidemiologists for the analysis of co- 
hort data are conditional, association models (e.g., logistic, pooled logistic. Cox 
proportional hazards). The time-dependence of exposure effects usually are ad- 
dressed in one of two ways: 1) ignored, in that only baseline exposure is consid- 
ered or the time-dependent exposures are summarized in one cumulative mea- 
sure whose effect is assumed to be confounded by baseline variables only (e.g. 
"pack-years" is a cumulative measure of exposure to cigarette smoke commonly 
used in occupational epidemiology) ; and 2) risk-set sampling, in which the ex- 
posure measure (possibly cumulative as just described) is updated at specific 
time points and exposure effects estimated based on the probability of outcome 
within groups classified by exposure in the current or most recent time interval. 
Thus, the full complexity of the effect of an exposure history on a specific health 
outcome is typically lost in conditional, association models. More importantly, 
these models also suffer from a major drawback: they typically produce biased 
effect estimates in scenarios involving time-dependent confounders also "on one 
of the causal pathways" of interest (3, 13) ■ Finally, these models do not provide 
direct population-level estimates of exposure effects, which often are of most rel- 
evance to public health. Only conditional associations between the exposure and 
outcome can be directly derived from traditional regression models. Indirectly 
however, causal inferences can be drawn from such models in particular cases 
and under additional assumptions. In such scenarios, it is important to note 
that only conditional (i.e. not population-level) effect estimates can typically be 
derived from such models. 

MSMs define parameters with a direct population-level causal interpreta- 
tion. Moreover, estimators of MSM parameters can account properly for time- 
dependent confounders also on one of the causal pathways of interest between 
the exposure and health outcome (fioh. Thus, these models can represent accu- 
rately population-level effects of histories of exposures on the health outcomes: 
the cumulating effects of exposures as well as acute effects from exposures expe- 
rienced over shorter time intervals. This is particularly important for exposures, 
such as ambient air pollutants, whose acute and chronic exposure effects lead 
to different health outcomes and/or contribute together in the occurrence of 
serious health outcomes such as heart attacks and death due to diseases of the 
heart and lung (|l|). 
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The current MSM methodology represents exposure effects by considering 
the exposure histories over the entire follow-up intervals that precede the oc- 
currence of time-specific health outcomes. In cases where consideration of the 
entire time interval prior to an event occurrence may not be relevant, based on 
subject matter, this omnibus handling of time possess an important limitation. 
This paper presents an extension of current MSM methodology to allow for 
more flexible analysis of time effects of exposure, based on a priori or ad hoc 
considerations of specific periods of time antecedent to an event. 

Note that this modeling extension was originally proposed and motivated 
in earlier work with discretized failure-time data and called partially marginal 
structural models (jH, S)- We develop here the formal statistical framework be- 
hind these models, compare it to the MSM framework and further motivate this 
approach with two other research problems. 

First, we undertook a study to determine the extent to which reductions 
in ambient air pollution consequent to regulations propagated since 1980 by 
the California Air Resources Board to reduce air pollution in the Los Angeles 
(LA) Basin lead to measurable health benefits. The basic time unit for the data 
was the quarter (3 months), and we had 84 such time units. The geographic 
area of interest was divided into 195 10 x 10 km grids, based on know patterns 
of air pollutants and meteorology in the LA Basin. The exposures of interest 
were quarterly concentrations of ambient air pollutants (e.g., ozone, oxides of 
nitrogen, particulate matter with a mass medium aerodynamic diameter of 10 
microns or less). A variety of health outcomes are being considered: quarterly 
hospital discharges and mortality rates for various chronic lung and cardiovas- 
cular diseases. Population denominators and demographic data are available 
on a quarter-spatial unit-specific basis. Over the 20 years encompassed by the 
study, there have been large temporal trends in demography that have lead to 
changes in population susceptibility to certain diseases of interest. For example, 
there has been a large infiux of Mexicans into the study area. Mexicans are 
known to have decreased risks for asthma and increased risks for diabetes mel- 
litus, an important underlying risk factor for cardiovascular disease. Moreover, 
since many of these immigrants are of low socioeconomic status, they may be 
more likely to live in closer proximity to sources of ambient pollutants (i.e., near 
Freeways). Thus, demography is an important temporal confounder. Changes 
in medical care over the study period also has affected the occurrence of health 
outcomes. In our descriptive analyses, important temporal trends for hospital 
discharges for asthma, chronic obstructive lung disease and various cardiovas- 
cular disease were observed. Temporal trends of disease-specific mortality are 
expected, e.g. there has been a decline in age-specific death rates from specific 
heart disease due to improvements in medical care. Based on the above, "time" 
becomes an important confounding variable to capture all of the unmeasured 
temporally-related factors that we have not measured and the residual temporal 
confounding for those factors that we have measured. Since our initial analysis 
focused on hospital discharges for asthma in children ages birth to 19 years, a 
central issue that emerged was the relevant exposure time for investigation of 
the exposure effect on the outcome rates. Based on available data on the ef- 
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fects of ambient air pollution on hospital discharges for asthma, it did not seem 
reasonable to extend the exposure period much beyond the 12 months prior to 
a given quarter. Consequently, there was a need to modify the how exposure 
histories are handled in the existing MSM methodology. 

Second, the proposed HRMSM methodology also has application to panel 
data that are being collected as part of a study of the relation between responses 
to short-term increases in ambient air pollutants and the long-term changes 
in symptoms and disease severity in children with asthma (Fresno Asthmatic 
Children's Environment Study - FACES). In this study, subjects participate in 
multiple 14-day panels during which time each subjects provides daily data on 
lung function, respiratory symptoms and daily activities. Analyses focus on the 
causal relation between daily symptoms or lung function and exposure to one 
or another pollutant over one or more days prior outcome report. Confounders 
for these analyses relate to meteorology and the effects of other pollutants not 
of primary interest in a given analysis. Virtually all studies to date evaluate the 
temporal relation between pollutants and symptoms/lung function through one 
of several modeling techniques based on traditional regression models where the 
study of the exposure effect is limited to: 1) specific lag effects (e.g., exposure on 
one or more days prior to outcome report) ; 2) specific average-over-several-days 
effects, or 3) polynomial-distributed-lag- functions effects (0). These approaches 
underscore the implicit understanding that the exposure effects of interest are 
limited in time. 

To reflect these subject-matter considerations, we propose HRMSMs as an 
extension of MSMs to provide greater flexibility in the evaluation of exposure 
effects over time within a rigorous statistical framework for causal inference 
based on counterfactuals. Once the investigator specifies the time frame over 
which the effect of the pollutant exposure is of interest, a recently developed 
data-adaptive methodology for model selection (jiol ) can provide guidance on 
how each level of the pollutant during the specified time frame should enter the 
HRMSM (e.g. moving averages or lags) that represents the hypothesized causal 
relationship between the pollutant and asthma outcome. 

In section[2l we introduce HRMSMs and in section [3] present three estimators 
of HRMSM parameters developed based on a so called i-specific counterfactual 
framework. HRMSMs are then compared to MSMs in section 2] before illus- 
trating the application of the proposed HRMSM methodology with a real-life 
application in section [5] Finally, we summarize and discuss the results presented 
in this article in section [6l For completeness, we include technical details in a 
comprehensive appendix. 

2. History-Restricted Marginal Structural Model 
2. 1 . Data structure and question of interest 

For all experimental units in a random population sample of size n, we observe 
a treatment regimen (^(0), . . . , A{K)) over time t = 0, . . . , if and a covariate 
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process (i(0), . . . , L{K +1)) measured at baseline and after a new treatment is 
assigned. The covariate L{t) is measured after A{t — 1) and before A{t). Note 
that K -\- \ represents the length of the treatment regimen in the appropriate 
unit of time and n the sample size. 

In the formal counterfactual framework for longitudinal studies (26h. the data 
are represented as n independent and identically distributed (i.i.d) realizations 
of: 

O = (L(0), ^(0), L(l), A(l), . . . , L(K),A(K),L{K+l)) = {A{K),L{K+1)) ^ P, 

where P represents the distribution of the stochastic process O, referred to 
as the observed data, and the general notation ~{t) represents the history of 
the variable '•' between time and t: a) ~{t) — (-(0), . . . , -{t)) if t > and b) 
~{t) is empty \i t < 0. We extend this notation with the notation ~(t_,t+) to 
represent the history of the variable ■ between time points t_ and i-|_: where a) 
~(t_, i+) — (-(t-), . • . , if t- < i+, and b) ~(t-, i+) is empty otherwise. We 

thus have"(t) =~{0,t). 

We denote with s an integer between and K + I specified by the inves- 
tigators to represent the history size for which the exposure effect is of in- 
terest in the analysis. For s ~ 1 < t < K, we define V{t — s + 1) as a sub- 
set of {A{t — s),L{t — s + 1)) to represent time-specific baseline variables, i.e. 
V{t — s -I- 1) C {A{t — s), L{t — s + 1)) and V{t — s + 1) represents the baseline 
variables if i — s -I- 1 is considered the study's starting time, i.e. the time from 
when the treatments and outcomes are considered of interest. In particular, we 
define a subset of the 'true' baseline variables as = V^(0) C L{0). We de- 
note the time-dependent outcome with Y{t): Y{t) e L{t), and T represents the 
set of time points t such that the outcome, Y{t + 1), is of interest. We have 
T C {0, . . . , K}. Typically T — {0, . . . , K} except when one is interested in the 
outcome collected at the end of the study only, i.e. when T = {K}. 

The question of interest is to investigate the causal effect of the treatment on 
the time-dependent outcomes of interest. In the literature, this problem has been 
addressed with MSMs. We propose in this article to address the same problem 
with the proposed HRMSMs. In MSM-based causal inference, the investigation 
of the causal relationship of interest relies on a representation of the effects of 
the treatment history between baseline and time point t, A{t), on the time- 
dependent outcome, Y{t + 1), for all i G T (see figure p. In HRMSM-based 
causal inference however, the investigation of the causal relationship of interest 
relies on a representation of different causal effects: the effects of the treatment 
history between time points t — s + 1 and t, A{t — s + l,t), on the time-dependent 
outcome, Y{t+ 1), for all t ^ T. Compared to MSM-based causal inference, the 
effect of the treatment is thus investigated for a history of treatment that is 
restricted by the investigators based on considerations discussed later in this 
article (see figure [2|). In other words, MSMs and HRMSMs can be viewed as 
two different statistical alternatives for the investigation of any given causal 
relationship, each providing different information and description of the causal 
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effect of interest. We argue that an HRMSM-based causal inference strategy 
may often be more suitable or at least more practicable than an MSM-based 
causal inference strategy for public health research. 

2.2. Assumptions 

Existence of counterfactuals: We assume the existence of the following 
treatment-specific processes, also referred to as a counterfactual processes 
[20I ). La(K){K + 1) for every treatment regimen a{K) = [a(0) , . . . , a(K)) G 
Av (K) , where Av (K) designates all possible treatment regimens between time 
points and K as a function of the baseline variable V only. In other words, 
Av{K) is the support of the conditional distribution of A{K) given V, g{A{K) \ 
V). We denote the so-called full data process with X and and its distribution 

with Fx: X = {V, + Fx. 

Note that the existence of the counterfactual process La(K){K + 1) for every 
treatment regimen a{K) G Av{K) implies the existence of the counterfactual 
processes La(t){t + 1) = ia(t),A(t+i,A:) + I) C X for every t = 0,...,K -1 
and every treatment regimen a{t) — (a(0), . . . , a{t)) G Av{t) where Av{t) des- 
ignates all possible treatment regimens between time points and t, i.e. the 
support of the conditional distribution of A{t) given V, g{A{t) \ V). We have 
Av{t) = {a{t) : 3 a'{K) G Av{K) a{t) = a {t)} for t = 0, . . . , K - 1 and Av{t) 
is thus entirely defined by Av{K)- 

Consistency assumption: At any time point t, we assume the following link 
between the observed data and the counterfactuals: L{t) = L^(^j^-^{t). Under this 
assumption, we have: O = {A{K),L^(^j^^{K + 1)) = (f){A{K),X), where ^ is a 
specified function of the full data process X. This notation indicates that the 
problem can be treated as a missing data problem. Only the treatment-specific 
process associated with the observed treatment A{K) is observed; the others 
are missing. 

Temporal Ordering assumption: At any time point we assume that any 
treatment-specific variable can only be affected by past treatments: = 
La(t-i){t) for t — 0, . . . , if -I- 1, where La(_i)(0) = L{Q). This assumption is 
typically implied by the data collection procedure: the covariate L{t) (which 
contains the outcome at time t) is measured after A{t — 1) and before A{t). 

Sequential Randomization Assumption (SRA): At any time point t, we 
assume that the observed treatment is independent of the full data, given the 
data observed up to time point t: A{t) _L X \ A{t — 1), L{t). Under the SRA, the 
treatment mechanism, i.e. the conditional density or probability of A{K) given 
X: g{A{K) \ X), becomes: 



K K 

g{A{K) I X) = X{g{A{t) \ A{t -l).X) l[g{A{t) \ A{t - l),L{t)). 
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The SRA implies coarsening at random ([3) and, thus, the likehhood of the 
observed data factorizes into two parts: a so-caUed Fx and g part. The Fx 
part of the hkehhood only depends on the full data process distribution, and 
the g part of the likelihood only depends on the treatment mechanism. As 
a consequence of this factorization of the likelihood under the SRA, we now 
denote the distribution of the observed data with Pf^.s ^^i^ the likelihood of O 
is: 

Fx part 

' ^ q part 

K+l 

C{0) /(L(0)) n fim I Ut 1), ^(^ - 1)) 9{A{K) I X) . 

Qfx 

In addition, we denote the set of conditional densities or probabilities that define 
the Fx part of the likehhood, except for f{L{0)) with Qfx- 

2.3. HRMSM and causal parameter of interest 

We define an HRMSM as a model for a feature (e.g. expectation) of the distri- 
bution of the counterfactual outcomes, i^A(t-s).a(t-s+i.t)(^ + 1)' corresponding 
with treatment interventions between time point t — s + 1 and t only (treat- 
ments between time points and t — s, A{t — s), are left random), possibly 
conditional on the baseline covariates at time t — s + 1, V{t — s + 1), for all pos- 
sible treatment interventions d{t — s + l,t) and t ^ where Tg represents the 
set of time points t such that the outcome y(t + 1) is of interest and < > s — 1, 
%^ {t : t > s - 1}. We have C {s - 1, . . . , K]. Typically we will have 
% = {s-l,...,K}. 

In section 21 we discuss the interpretation of HRMSM parameters and how 
they represent the causal relationship of interest for a given value for s compared 
to MSM parameters. By convention in this article, the random portion of the 
treatment history defining counterfactuals is excluded from the counterfactual 
notation and thus we adopt the following notations: y^(t_5) a(t-s+i t)(^ + 1) = 
Fa(t_,+i,t)(t-hl). 

Typically and specifically in this article, one is interested in average causal 
effects per stratum V{t — s + 1) of the population which can be represented 
by causal parameters defined by HRMSMs for Epx ,g{Ya.(t-s+i,t){t + 1) I ^(^ — 
s + 1)) for t ^ Tg. Similar to MSM analyses, investigators can choose to model 
Epx ,g0^a{t-s+i,t){t + 1) I y{t ~ s + 1)) for each t G separately with a strat- 
ified model or simultaneously with a pooled model ([sl). We denote the causal 
parameters defined by a stratified and pooled HRMSM with Pt{Fx,g \ ■) (one 
parameter for each t) and l3{Fx,g \ ■) (single parameter for all t), respectively, 
to indicate that they are mappings from the space of distributions {Fx , g) to 
the space of real numbers and that these mappings are functions of modelling 
assumptions represented by Similar to the two MSM approaches developed 
in the literature , a parametric or nonparametric HRMSM approach can 

be adopted, each relying on different modelling assumptions. 
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Note that unlike the class of MSMs, HRMSMs are introduced as a class of 
mixed full and observed data models, since an HRMSM models the marginal 
distribution of counterfactuals, where part of the treatment history is left ran- 
dom. The distribution of the random portion of the treatment history is defined 
by the treatment mechanism, g, and, thus, is identified by the observed data. 

3. HRMSM estimation 

For clarity, we only present in this section a summary of the theoretical results 
detailed in the appendix as they relate to HRMSM estimation. For conciseness, 
we illustrate their practical applications with the example of one estimator only. 

3.1. The t-specific counter factual framework 

The time-specific (t-specific) counterfactual framework can be viewed as an 
extension of the conventional counterfactual framework on which MSM-based 
causal inference relies. This latter mathematical construct was described in sec- 
tions 12.11 and 12.21 It provides the rigorous framework to define, identify and 
estimate MSM parameters with the full and observed data based on a sufficient 
set of assumptions presented in section [2T2] We describe the t-specific counter- 
factual framework in detail in appendix [A] and introduce it because it allows us 
to generalize the MSM estimation procedures to HRMSM estimation procedures 
with minimum effort. After linking the t-specific counterfactual framework to 
the conventional counterfactual framework in appendix |X1 we show in appendix 
[B]that the sufficient set of assumptions for MSM estimation presented in section 
12.21 is also sufficient for HRMSM estimation. An important consequence is that 
the choice of HRMSM-based causal analysis over MSM-based causal analysis 
becomes only a matter of practical considerations (statistical power and com- 
puting issues) and, above all, subject-matter considerations (the relevance of 
the causal effect representation to public health research). These considerations 
are developed in section [H 

3.2. Link between the conventional and t-specific counterfactual 
frameworks 

Figure [5] in the appendix illustrates based on an example of a longitudinal 
study with short follow-up the link between the longitudinal data representa- 
tion in the conventional counterfactual framework and its representation in the 
time-specific counterfactual framework. Note that in the conventional counter- 
factual framework the data are approached as a single entity, O, in the sense 
that the treatment is defined once and for all as a history A{K) and the outcome 
is time-dependent, Y{t) £ L{t). On the other hand, in the t-specific counter- 
factual framework the data are viewed as layers of separate entities, O*. For 
each O*, the treatment and outcome of interest are redefined along with the 
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baseline covariates (highlighted in yellow/light gray on figure [5|). Note that for 
each O* , the outcome is no longer time-dependent but correspond with the last 
outcome collected, F* £ + 1), (highlighted in orange/dark gray on figure 
[5]) and the treatment history size is fixed to a user-specified value s = 2. In the 
conventional counterfactual framework, the investigator examines the effect of 
A{K) on Y{t) for alH G T based on MSMs for the full data associated with O 
whereas in the ^-specific counterfactual framework, the investigator can examine 
the effect oi A\t - s + l,t) on F* for t e % based on MSMs for the t-specific 
full data associated with O* . Figure [5] illustrates how the ^-specific counter- 
factual framework can be viewed as a collection of conventional counterfactual 
sub-frameworks with distinct definition of the outcome, treatment and baseline 
covariates. These conventional counterfactual sub-frameworks differ from the 
conventional counterfactual framework in the sense that the treatment history 
is of size s ^ K + \ and the outcome of interest is no longer time-dependent. 



3. 3. Three estimators of HRMSM parameters 

We show in the appendix that an HRMSM corresponds to a t-specific MSM and 
the intersection of i-specific MSMs in a stratified and pooled analysis, respec- 
tively. This result implies that estimating functions for HRMSM parameters are 
t-specific MSM estimating functions (stratified analysis) or sums of t-specific 
MSM estimating functions (pooled analysis). Thus, we can easily extend the 
three estimators developed for MSM parameters to HRMSM parameters: the 
Inverse Probability of Treatment Weighted 17 Jl ^ . the G-computation 



dil [H, i, [H, [H, lii, i, i) and Double Robustllinifllj estimators. To illus- 
trate how the definition and implementation of MSM estimators are extended to 
HRMSM parameters based on the theoretical results in the appendix, we focus 
on IPTW estimation which has been favored in real-life applications. 



3.4- IPTW estimator: definition and implementation 

We present the two classes of IPTW estimators for HRMSM parameters in the 
stratified and pooled analyses based on a parametric HRMSM approach. The 
reader can easily extend these estimators to the nonparametric approach based 
on the MSM literature (0,|1). 



3.4- 1- Stratified analysis 

In this analysis, causal effects are modelled separately for each time point t Ts, 
i.e., one separately investigates the causal effects on the outcomes of interest, 
Y{t+1) ioT t gTs, through the estimation of distinct causal parameters Pt{Fx \ 
■) ioT t gTs defined based on I = Card(7^) distinct models rrit for < G 7^. 

Under this model, the investigation of the causal effects of interest based 
on parametric models corresponds to the estimation of the / distinct causal 
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parameters (31 = Pi{Fx,g \ ™t,5) defined such that: 

^Fx,9(>"a(i-s+i,t)(t+l) I V{t-s+l)) - mt{a{t-s+l,t),V{t-s+l) \ /3*) for t E %. 

For a given t, the ciass of IPTW estimating functions for (31 with nuisance 
parameter g is defined as: 

DiAu \ n. 3) — — ; = : n Gn>, 

uU-s+i9iMj)\Au-i),L{j)) y 

where: 

• e(^) = Y(t + 1) - mt(i(i - s + 1, i), y(t - s + 1) | /3), 

• 7i is a set of non-null functions of A(t — s + 1, t) and V(t — s + 1). 

We denote estimators of g and /i with g^ and /i„. 

We define the Experimental Treatment Assignment (ETA) assumption (0) 
for estimating /9j* as follows: 

max — — — — — < oo i^x - a.e, 

a(t- X)^Av(t- ,t) * 

J] g(a(j) |A(t-s),a(t',j-l),i(j)) 

where t~ =t — s + 1 and Ay (i — s + 1 . t) is the set of possible treatment regimens 
between time points i— s+l and t, i.e. the support of the conditional distribution 
of A(t - s + 1, t) given A{t - s) and V, g(^A(t - s ^ \,t) \ A(t - s),V). 

The IPTW estimator of (31 is defined as the solution of the estimating equa- 
tion associated with the observed data O and the IPTW estimating function at 

n 

Y^^hM l5n,/3) = 
1=1 

Under regularity conditions, the IPTW estimator of (31 is consistent and asymp- 
totically linear if the ETA assumption holds and if gn is a consistent estimator 
of g. 

In practice, the implementation of the IPTW estimator of the HRMSM pa- 
rameter (31 is identical to the implementation of MSM parameters in a stratified 
analysis with the difference that the treatment of interest for outcome Y(t^V) is 
limited to A{t — s -f 1, t) instead of A{t). In other words, the IPTW estimate can 
be obtained with a weighted least squares regression oiY(t + \) on A(t — s + \^ t) 
and V{t — s + 1) based on the parametric model rrit and weights inversely pro- 
portional to the estimated treatment mechanism. Like in MSM estimation, the 
numerator of the weights is defined by the choice for h (e.g. stabilized versus 
unstabilized weights (|l9r )). Unlike MSM estimation, note that the treatment 
mechanism is limited to treatments assigned between time points t — s + 1 and 
t: A{t — s + l,t) and not all time points between baseline and t: A{t). 
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3.4-2. Pooled analysis 

In this analysis, causal effects are modelled simultaneously for each time point 
t G Ts, i.e. the change of the causal effect on the outcome over time is represented 
by a smooth function of time: one simultaneously investigates the causal effects 
on Y{t + 1) for t G 7^, through the estimation of a single causal parameter 
f3{Fx,g I •) defined based on a single model m{t, a{t — s + 1, t), V{t — s + 1) \ f3). 

Under this model, the investigation of the causal effects of interest based on 
parametric models corresponds to the estimation of the single causal parameter 
P* = p*{Fx,g I rn,g) defined such that 

Ep^jY-^it-s+iMt + 1) I yit-'^ + 1)) - ™(^, a{t-s + 1, t),V{t - s + 1) I /3*) 
for t € %. 

For a given t, the class of IPTW estimating functions for /3* with nuisance 
parameter g is defined as: 

ir,^t-s+i9iMj)\AU~l),Lij)) 
where: 

• £{/3) = Y{t + 1) - m(t, A{t - s + l,t),V{t -s + l)\ f3), 

• 7i is a set of non-null functions of t, A{t — ,s + 1, t) and V{t — s + 1). 

We denote estimators of g and h with gn and /i„. 

We define the Experimental Treatment Assignment (ETA) assumption (0) 
for estimating /3* as follows: 

h{t,d{t-,t),V{t-)) 
max max < oo rx — a-e, 

t£T, a{t-,t)GAv(t-,t) * _ 

n g{a{j)\A{t~s),a{t-,j-l),L{j)) 

3=t- 

where for each t, t~ = t — s + 1 and Av{t — s + l, t) is the set of possible treatment 
regimens between time points t — s + \ and t, i.e. the support of the conditional 
distribution oi A{t-s + l,t) given A(t-s) and V , g{A{t- s + l,t) \ A{t-s),V). 

The IPTW estimator of /3* is defined as the solution of the estimating equa- 
tion associated with the observed data O and the IPTW estimating function at 
gn ■■ 



Under regularity conditions, the IPTW estimator of (3* is consistent and asymp- 
totically linear if the ETA assumption holds and if g„ is a consistent estimator 
of g. 

In practice, the implementation of the IPTW estimator of the HRMSM pa- 
rameter P* is identical to the implementation of MSM parameters in a pooled 
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analysis with the difference that the treatment of interest for outcome Y{t + 1) 
is hmited to A{t — s + 1, t) instead of A{t). In other words, the IPTW estimate 
can be obtained with a pooled weighted least squares regression of Y{t + 1) for 
alH S 7^ on A{t — s + 1, t) and V{t — s + 1) based on the parametric model m 
and weights inversely proportional to the estimated treatment mechanism. Like 
in MSM estimation, the numerator of the weights is defined by the choice for h 
(e.g. stabilized versus unstabilized weights (jlOj)). Unlike MSM estimation, note 
that the treatment mechanism is limited to treatments assigned between time 
points t — s + 1 and t: A{t — s + l,t) and not the entire treatments between 
baseline and time point t: A{t — s). 



4. When and why prefer HRMSM-based versus MSM-based causal 
inference in practice? 

4.1. MSM parameters: interpretation and causal effect 
representation 



MSMs were introduced as a class of full data models which define parameters 
based on a feature of the marginal distribution of the following counterfactual 
outcomes: Ya(t){t+1) possibly conditional on the baseline covariates V . Typically 
and specifically in this article, one is interested in average causal effects per 
stratum V of the population which can be represented by causal parameters 
defined by MSMs for EFx{Ya(^t){t + I) \ V) for t e %. We denote a causal 
parameter defined by an MSM with (3t{Fx I ■) to indicate that it is a mapping 
from the space of full data distribution Fx to the space of real numbers and 
that this mapping is a function of modeling assumptions represented by •. 

Two approaches to causal inference based on MSM have been proposed. They 
provide different representations of causal effects with distinct causal parame- 
ters. Initially, a parametric MSM approach to causal inference was developed 
(jl4l ) that relies on correct specification of a parametric MSM. Recently, a new 
approach based on nonparametric MSM was introduced that does not 

require to assume a correctly specified MSM and that generalizes the definition 
of causal parameters. This later approach is more realistic if one believes that 
correct specification of a parametric MSM is unlikely in practice. 

In addition, both MSM approaches can be based on either a stratified or 
a pooled analysis, i.e. distinct models, 'mt{a{t) ,V \ (it) , or a single model, 
m{t,a{t),V\ /3), for EF^{Y-,^^t){t + 1) | for t G T i). 

Independently of the MSM approach chosen (nonparametric versus paramet- 
ric and pooled versus stratified) , MSM parameters represent the causal effects of 
the treatment histories, A{t), on the outcomes, Y{t+1), for t G T. Note that this 
implies that in MSM-based causal inference, the causal effect of the treatment 
on the outcome collected at time point t is always investigated for a treatment 
history of size t. As a result, the causal effect on the outcome collected at time 
t, Y{t), is defined based on larger treatment histories as t increases, i.e. as the 
outcome is collected later in the longitudinal study. This feature of this causal 
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Fig 1: Illustration of the MSM representation of causal effects in a longitudinal 
study with time-dependent outcomes. Each ellipse indicates the relevant time- 
span over which the treatment's effect on a given outcome is investigated. Each 
cross indicates an outcome and each arrow represents one of the studied effects 
between treatment histories (ellipses) and outcomes (crosses). 



analysis is illustrated in figure [T] The figure illustrates how causal effects are 
investigated in practice based on an MSM approach for a study where K — b, 
i.e. where the observed data is: 

O = (L(0), A(0), L(l), ^(1), L(2), A(2), L(3), A(3), i(4), A(4), L(5), A(5), L(6)). 

For instance, the causal effect of the treatment on the outcome collected at time 
point f = 6 is investigated for a treatment history of size 6, ^(5), as represented 
on figure [T] by the ellipse that covers time point to 5 and the arrow that 
connects the ellipse to the outcome collected at time point 6 represented by a 
cross. 



The investigation of causal effects with MSMs (see figure[T]) raises three potential 
concerns that are likely to become more significant for longitudinal studies with 
longer follow-ups: 

• computational intractability when proceeding to MSM estimation 

• a disbelief about the subject-matter relevance of the causal effects inves- 
tigated 

• a statistical power problem 

The first issue is best illustrated with the implementation of the G-computa- 
tion estimator and was fully developed in previous work (0, [3). More generally, 
implementation of MSM estimators, like the IPTW estimator, becomes less 
practicable as the follow-up time increases. 

The second issue can easily be illustrated with the following hypothetical 
study. Consider a longitudinal study during which individuals are treated or not 
every day over three months (90 days) with a new medication for headache relief 
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and monitored for headache symptoms. Now consider the last outcome, y(90), 
collected after a treatment history of 89 days, A{89). In MSM-based causal 
inference, the investigation of treatment effect on the last outcome measured 
would be based on the estimation of a causal parameter representing the effect 
of the treatment history A(89) on F(90). Most investigators would argue that 
looking a the effect of such a medication taken 3 months prior outcome report 
is likely to be of little interest since: 1) the relief effect of such medication 
usually does not carry over such a long period of time; and 2) because the 
drug effect that is pursued, i.e. of interest, for such a treatment is a short- 
term relief. In other words, investigation of the effect of a headache reliever 
absorbed 3 months prior outcome report is not of primary interest. It may also 
often make no sense to study such a long lag effect in practice if the treatment 
is known to act over a short-term time scale exclusively. Note that a natural 
way to overcome this problem in practice would be to investigate short-term 
effects based on a parametric MSMs only involving the latest treatments received 
before the outcome is collected. For example a parametric MSM for ¥3(39) (90) 
may only rely on the last two treatments absorbed to explain the outcome, 
e.g. £;(ra(89)(90)) = po+ /3ia(88) + f32a{89). Note however that such an MSM 
assumes that treatment before time 88 has no effect on the treatment which 
may be incorrect and estimation of f3i and P2 should still rely on the complete 
treatment history between time point and 89. We argue in this article that 
the MSM approach can be improved to better identify causal effects that are 
truly of interest from a subject-matter point of view without making restrictive 
assumptions often incorrect in practice. The HRMSM proposed in this article 
addresses this issue and mitigates the other two concerns discussed here. 

The third issue is not illustrated in this article with a concrete example but 
has been discussed in previous work S) and is related to the Experimental 
Treatment Assignment assumption (jTI). In short, in statistical analyses based 
on MSMs, the longer the treatment history, i.e. the study follow-up, the more 
complex the description of the treatment effect becomes. Thus, it is likely that 
the information required to understand long-term effects will be very impor- 
tant and beyond the reach of most investigators. Concretely, the decrease in 
statistical power can be explained by an increasing practical violation of the 
ETA assumption as the treatment history increases. As a result, even when it 
is sensible to investigate causal effects based on MSMs in longitudinal studies 
with long follow-up, investigators may still wish to revise their study aims and 
lower their research ambitions for the sake of practicability. Note that it is this 
statistical power issue which initially motivated the introduction of HRMSMs 
with failure-time data (0, 0)- The authors indeed noted that MSM estimation 
with the IPTW estimator suffered from a lack of precision for long treatment 
histories due to an increasing variability in the weights. This phenomemon can 
be explained by practical violation of the ETA assumption when the treatment 
history is long (|7|). 
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4-2. HRMSM parameters: interpretation and causal effect 
representation 

Both the parametric and nonparametric MSM (14, 1^, d) approaches that have 



been proposed in causal inference can be directly extended to HRMSM-based ap- 
proaches. Similarly, the corresponding parametric and nonparametric HRMSM 
approaches to causal inference provide different representations of causal effects 
with distinct causal parameters. 

In addition, both HRMSM approaches can be based on either a stratified or a 
pooled analysis, i.e. distinct models, mt{a{t~s+l, t), V{t — s+l) \ (3t), or a single 
model, mit ait - s + I, t),V{t ~ s + 1) \ P)Jor Ep^ ,g{Ya(t^s+i.t) I V{t - s + 1)) 

for te%m- 

Independently of the HRMSM approach chosen (nonparametric versus para- 
metric and pooled versus stratified), the causal effect of the treatment on the 
outcome is always investigated for a treatment history of fixed size, s, in HRMSM- 
based causal inference. As a result, the causal effect on the outcome collected at 
time t, Y{t)^ is defined based on a fixed treatment history even as t increases, 
i.e. as the outcome is collected later in the longitudinal study. This feature of 
this causal analysis is illustrated in figure [2l It illustrates how causal effects are 
investigated in practice with HRMSMs for a study where K = b and s = 2, i.e. 
where the data collected is: 

(i(0), A(0), L(l), A(l), i(2), A(2), L(3), A(3), i(4), A(4), L(5), A(5), i(6)). 

For instance, the causal effect of the treatment on the outcome collected at 
time point i = 6 is investigated for a treatment history of size 2: (A(4), A(5)) as 
represented on figure[2]by the ellipse that covers time points 4 to 5 and the arrow 
that connects the ellipse to the outcome collected at time point 6 represented 
by a cross. 



The following three practical considerations may often lead investigators to pre- 
fer a causal analysis based on HRMSMs above an analysis based on MSMs for 
causal inference problems with longitudinal data: 

• Computational tractability, 

• Causal effect representation that is most relevant to public health research, 

• Statistical power. 

Support for these claims can be found in section IT] and 14 . 1 1 where the limitations 
of MSM-based causal inference are underscored. 

5. Illustration with an air pollution study 
5.1. Specific aims and data 



We illustrate the application of HRMSMs with one of the air pollution study 
already described in section [TJ Its primary goal is to investigate the extent to 



R. Neugebauer et al. /Causal inference with HRMSMs 



134 



t=0 12 3 4 5 6 







/ 
I 

V 






\ /N 
\ y 

< X 

/A \ , 











Fig 2: Illustration of the HRMSM representation of causal effects in a longi- 
tudinal study with time-dependent outcomes and s — 2. Each ellipse indicates 
the relevant time-span over which the treatment's effect on a given outcome is 
investigated. Each cross indicates an outcome and each arrow represents one of 
the studied effects between treatment histories (ellipses) and outcomes (crosses) . 



which reductions in ambient air pollution consequent to regulations propagated 
since 1980 by the California Air Resources Board to reduce air pollution in the 
Los Angeles (LA) Basin lead to measurable health benefits. In this example, we 
will specifically focus on the reduction of ozone levels and its effect on health 
as measured by asthma-related hospital discharges in children ages birth to 19 
year^. 

The experimental units are 10 x 10 km geographical areas of the LA Basin. 
Data were assembled on a total of n = 195 such grids over 72 [K ~ 70) quar- 
ters to answer the aforementioned question of interest. For each grid, i, and 
each quarter t, the data set contains measurements for 1) the exposure of inter- 
est: ozone levels denoted with A(t), 2) the outcome of interest: asthma-related 
hospital discharge counts denoted with C(t) and the total number of hospital 
discharges N(t) in children ages birth to 19 years, and 3) 56 socioeconomic and 
demographic covariates denoted with L(i) (e.g. age, income, racial and gender 
structures as measured by ratios of the grid census). 

5.2. Model and parameter of interest 

To address the question of interest, we believe that an MSM-based statistical 
strategy to analyze these data is most natural. Indeed, the overall subject-matter 
objective is two-fold: 1) we wish to verify that exposure to ozone indeed lead to 
averse public health impact, i.e. we wish to evaluate the causal effect of ozone 
on the asthma-related hospital discharge rate; and 2) after verification of the 
averse nature of this effect, we wish to evaluate the public health benefits related 
to (caused by) decrease in ozone levels over 20 years (see figure [3]), i.e. we wish 



^Results are for heuristic purposes only and are not presented as a definitive analysis on 
the subject. 
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Fig 3: Plot representing the decrease of ozone levels over 72 quarters for the 195 
grid units. The linear curve corresponds with the fit of a linear regression model 
of ozone levels against time. 



to contrast what would the population outcomes be had the ozone levels not 
decreased over time. 

A central issue that emerges is the relevant exposure time for investigation of 
the exposure effect on the outcome rate: it does not seem reasonable to extend 
the exposure period much beyond the 12 months prior to a given quarter. We 
indeed believe that most of the effect of ozone can be captured within this time 
frame. This choice of a time frame is reinforced by the fact that the experimental 
units are not individuals but geographical areas: the population constituting 
each grid is constantly changing and major modifications in the population 
structures occurred over the 20-year period studied. Our focus is, thus, on the 
sub-acute effects: we consider the effect of one-year exposure to ozone during 
which time the population in each grid is relatively constant. Even, if the effects 
of ozone were to extend much beyond this time frame, we believe that we would 
lack power to identify such effects with our finite data (n = 195). In addition, 
software development is too limited to date to engage in a more complex analysis 
where we would investigate the ozone levels over large histories. Consequently, 
an HRMSM is an appealing statistical tool for the analyses of these data. 

We propose to use a pooled, binomial logistic HRMSM with s = 4 (i.e., we 
consider the effect of 4 quarters of exposure to ozone) to represent the effect 
of interest in this analysis, i.e. we assume that the conditional distribution of 
Ca{t-3,t)it + 1) conditional on -/Va(t-3,t) (i + 1) follows a binomial distribution: 
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Ca(t-3,t)(i + l) ^S(7Va(t_3,t)(t + l),pa(4_3,4)(t + l)) with: 

/ g.(.-3,o(t + l) I ( A 

V^a(t-3,t)(i + l) ' ^ ^ 7 

= logit(^/3o + /?iniean(a(t - 3, t)) + /32/i(t) + /33/2(i) 

+/34mean(a(i — 3, + /35niean(a(t — 3,t))/2(i) 

+/36mean(a(i - 3, i))/i(i)/2(t) 



where: 



t e = {3, . . . , 70} 
logit(x) 



l+cxp (-(k)) ' 

. mean(a(t - 3, t)) = a(«-3)+a(t-2)+a(.-i)+a(t) ^ 

• /i(t) is a mapping from i to the year number associated with quarter i+ 1, 
and 

• /2 {t) is a mapping from t to the season associated with quarter t-\-\ (two 
6- month seasons are considered in the analysis only: winter and summer) . 

This model implies the following: 

log f T^^^^^^^^^^T^^) = /3o + /3imean(a(t - 3, t)) + Miit) + /33/2(i) 

+/34mean(a(i — 3, t))fi{t) + j3^mea,n{a{t — 3, t))f2{t) 
+Pe^ean{a{t - 3, i))/i(t)/2(i), 



because i? 



5.3. Estimation 

To fit the HRMSM, i.e. estimate (3 = (/3i, . . . , /Je), we chose the G-computation 
estimator. A detailed algorithm for the implementation of the G-computation 
estimator has been described in previous work (8, 9). This implementation re- 
lies, in practice, on a data reduction step based on a generalized linear model for 
the treatment mechanism. Based on the reduced data set, the G-computation 
estimate of (3 is obtained through Monte Carlo simulation based on a gener- 
alized linear model for the Qfx P^rt of the likelihood. All models involved in 
this estimation procedure were obtained through model selection. We utilized a 
recently developed data-adaptive model selection procedure (2^, [25l ) based on 
cross-validation and an aggressive model search algorithm known as the D/S/A 



algorithm (|2ll). 
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Fig 4: Effect of one-year exposure to ozone on astfima-related hospital discharges 
over time (curvihnear line). The straight lines represent the season-specific 
causal effects. Each point above the y = line is interpreted as a decrease 
in the odds of asthma-related hospital discharge incidence caused by a decrease 
of one unit in the average ozone quarterly levels over one year. 



5.4- Results and interpretation 

Note that the results presented are preliminarj0. Figure S] represents a sum- 
mary of the results with a plot of the temporal change of the effect of one-year 
exposure to ozone on the asthma-related hospital discharge rate as measured 
by $1 + P4fi{t) + ji^f^it) + P6fi{t)f2{t) where /3 represents the G-computation 
point estimate of f3. Each point above the y = line on figure |4] is interpreted 
as a decrease in the odds of asthma-related hospital discharge incidence caused 
by a decrease of one unit in the average ozone quarterly levels over one year. In 
addition, the unit effect of ozone on asthma-related hospital discharge incidence 
decreases linearly over time. Note also that a portion of the winter section of 
the curve is below the y = line which indicates a decrease in asthma-related 
admissions per unit increase of ozone levels. This decrease could represent un- 
controlled confounding, measurement errors and/or model misspecification. We 
did not explore this finding, since results have been presented to illustrate the 
method. Also, we do not provide confidence intervals which can be obtained 
with the bootstrap as is appropriate with the MSM methodology. 

^Results are for heuristic purposes only and are not presented as a definitive analysis on 
the subject. 
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6. Discussion 

In this paper, we have laid out the formal statistical framework for a new class 
of MSMs, HRMSMs. These models were initially introduced (0, [1) as alter- 
native causal inference tools to MSMs in order to improve statistical power 
when investigating the causal effects of a treatment history on a health out- 
come. We developed an extension of the conventional counterfactual framework 
that we called the t-specific counterfactual framework. This framework was in- 
troduced solely as a statistical artifice to provide the rigorous mathematical 
framework to develop consistent estimators of HRMSM parameters with mini- 
mal effort: the IPTW, G-computation and DR estimators. We have shown that 
these estimators of HRMSM parameters are consistent under the same model 
assumptions commonly adopted in the conventional counterfactual framework: 
existence of counterfactuals, consistency, time-ordering and sequential random- 
ization assumptions. 

In addition, we further argued in this article, based on practical considera- 
tions (computational tractability, relevance of the effect investigated, statistical 
power), that an HRMSM-based causal inference strategy may often be better 
suited for public health research than an MSM approach where the effect of 
an exposure history is modeled as a function of exposures experienced only in 
the last s periods prior outcome collection. Such an MSM approach implicitely 
relies on the assumption that exposures experienced before the last s periods 
had no effect on the outcome. This assumption typically does not hold in most 
applications. The HRMSM approach makes no such assumption, and so is more 
circumspect in its implications. We believe these considerations should motivate 
the application of this methodology in many epidemiological and clinical stud- 
ies. We now discuss the decision making about the history size, s, in real-life 
applications of the proposed HRMSM-based causal analysis. 

Decision about the value for the history size, s, should be based on the combi- 
nation of considerations about the analysis aims and a priori knowledge about 
the problem being studied as illustrated in section[5l This decision however can- 
not ignore practical considerations such as implementation issues and statistical 
power concerns. 

For instance, if a longitudinal study aims to investigate the causal effect of 
a new medication for headache relief whose action is likely not to carry over 
time beyond a few hours, then it will not make sense to choose a history size 
that extends well beyond the known lag effect of similar medication. Even when 
the treatment effect of interest is likely to carry over long periods of time, the 
subject-matter focus may be the investigation of short-term effects, in which 
case the investigators should consider small values for s. In addition, note that 
the larger s, the more complex the causal effect of interest that is captured in the 
time interval represented by s. As a result, the statistical power to investigate 
the causal effect of interest will likely decrease when considering larger history 
sizes s. Moreover, the larger s is, the less computationally tractable will HRMSM 
estimation be. 
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Nevertheless, choosing a suitable value for s based on these guidelines remain 
subjective and may lead to two scenarios: 1) the chosen history size, s, is larger 
than the maximum time interval over which the treatment of interest has an 
effect on the outcome; and 2) the chosen history size, s, is smaller than the 
maximum time interval over which the treatment of interest has an effect on 
the outcome. In the first case scenario, the model selection procedure for MSMs 
proposed by van der Laan and Dudoit (2003) ((23 ) can be used to identify the 
smaller component of the treatment history that is causally relevant. The second 
case scenario is most likely to occur in practice, since statistical power and im- 
plementation considerations often will prevent investigation of the causal effects 
of treatment histories that are too long. An HRMSM-based causal analysis still 
will provide valuable answers to the public health questions of interest, based 
on the available data even if the causal effect of the treatment over time will 
not be completely described (e.g. the maximum lag effect will remain unknown). 

We would also like to note that a variant of HRMSMs can be used for possibly 
discrete failure-time outcomes. For a discrete failure-time outcome, this would 
be a model for EFx,g{Ya_(^t-s+i,t){t + 1) 1 ^5(4-5+1,*) (0: V{t - s + 1)); this is a 
model for the probability of failure at time t + 1 under a given restricted treat- 
ment regimen conditional on having survived through time t under that regimen 
and possibly i-specific baseline covariates. Such models pose two problems for 
causal interpretation. First, as is generally true with comparisons of hazards, 
even in randomized trials, the comparison of two treatment regimens' effects 
at a given time point is not based on comparing outcomes for a common set 
of individuals; rather, it is based on comparing failure among survivors under 
one treatment regimen to failure for another set of survivors under the other 
treatment regimen. If subjects suceptible to harmful effects of a given treat- 
ment regimen fail early after being treated with that regimen, the depletion of 
susceptibles may lead to hazards under that treatment regimen being the same 
or even lower than in its absence. In simple randomized trials or with MSMs, 
this problem with interpretation can be mitigated by using the hazards under 
a given treatment regimen to compute and compare treatment-specific survival 
functions instead of treatment specific hazard functions. The comparisons of 
treatment-specific survival functions are comparisons for a common group of 
individuals under different treatment regimens, starting from the beginning of 
the regimen. In this variant of HRMSMs, treatment-specific survival curves can- 
not be computed, because unlike MSMs, hazards at different time points t in 
HRMSMs do not correspond with a single treatment-specific hazard of interest 
and can thus not be mapped into a treatment specific survival curve. This is 
particularly obvious if the HRMSMs condition on covariates V{t — s + 1) since 
the set of covariates V{t — s+1) will not be identical over time. These concern 
do not apply to HRMSMs for non-survival outcomes. 



Finally, we would like to clarify the nomenclature for models that have re- 
cently been proposed by different authors to essentially designate the same 
class of models: partially MSMs (0), History- Adjusted MSMs (jlj, and History- 



R. Neugebauer et al. /Causal inference with HRMSMs 



140 



Restricted MSMs. The difference in the nomenclature reflects the heterogeneity 
in the subject-matter problems for which different authors developed a model- 
ing approach based on the same class of MSMs. Unlike partially MSMs and 
HRMSMs, the primary motivation for the development of History-Adjusted 
MSMs was the inability to incorporate modification of the causal effect of treat- 
ment by time-varying covariates in MSMs for clinical decision making. It is in- 
teresting to note that solutions to these different problems can be found in the 
same class of MSMs. This confluence should further underscore the importance 
of the proposed modeling approach in practice. 

APPENDIX 

Appendix A: HRMSM estimation: the time-specific counterfactual 
framework 

In this section, we introduce the time-specific (i-specific) counterfactual frame- 
work which can be viewed as an extension of the conventional counterfactual 
framework on which MSM-based causal inference is based (see sections 12.11 and 
12. 2p . This latter mathematical construct provided the rigorous framework to 
define, identify and estimate MSM parameters with the full and observed data 
based on a sufficient set of assumptions developed in section [^21 We introduce 
the t-specific counterfactual framework because it allows us to generalize the 
MSM estimation procedures to HRMSM estimation procedures with minimum 
effort. In addition, we present in this section the sufficient assumptions for esti- 
mation of HRMSM parameters in the t-specific counterfactual framework. 

A.l. Data structures 

In this section, we adopt the notations introduced in the previous section to 
represent the treatments, covariates and outcomes collected at each time point 
t = 0, . . . , K+\ on each of the n experimental units: A{t), L{t), Y(t) respectively. 
We also adopt the notation V{t — s + 1) to designate a subset of {A{t — s),L{t — 
s + 1)). 

In the i-specific counterfactual framework, the representation of the data col- 
lected during a longitudinal study between time points and A' -I- 1 is based 
on a user-specified choice of a fixed treatment history size s > 0. We already 
discussed the interpretation of this parameter, s. In section |6l we discuss the 
decision making about its value in practice. 

In the t-specific counterfactual framework and for a given treatment history 
size s, the data are represented as n i.i.d realizations of K — s + 2 data structures: 

O* = {L\t- s + l),A\t- s + l),L\t- s + 2),A\t- s + 2),...,L\t), 

A\t),L\t + l)) 

= {A\t-s + 1, t),L\t - s + 1, i + 1)) - P*, 
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for all t such that s — 1 <t < K and where 

• P* represents the distribution of the stochastic process O* referred to as 

one of the t-specific observed data, 

• A*{t — s + l,t) represents the t-specific treatment process defined as 

A\j) = A{j) for aU j such that t- s + l<j <t, (1) 

• — s + 1, t + 1) represents the f-specific covariate process defined by: 

a) = for all j such that s + 1< j < i + 1 (2) 

b) L\t- s + l) = {A{Q,t- s),L{Q,t- s + 1)) (3) 
In other words, we have - s + 1, i + 1) = {A{t - s), L{t + 1)). 

Furthermore, we define V* as a subset of the baseline covariates in the f-specific 
observed data, O*, V* C L*{t — s + 1), such that: 

V* =V{t-s + l). (4) 

We define F* as the f-specific outcome of interest, F* G + 1), such that: 

Y* = Y{t+l). (5) 

In addition, we denote with Tg the set of time points t such that the out- 
come y* is of interest. We have Tg C {s — 1, . . . , K}. Typically we will have 
Tg = {s — 1, . . . ,K}. Note that the outcome of interest maybe time dependent, 
e.g. Tg = {s — 1, . . . , K}, while the single outcome of interest in the t-specific 
observed data is Y* = Y{t+ 1), i.e. Y{t) is not an outcome of interest in O* but 
only in 0*~^. 

Like in the conventional countcrfactual framework, the question of interest 
is to investigate the causal efi'ect of treatment A on the time-dependent out- 
come, Y G L. In the t-specific countcrfactual framework, this problem is ad- 
dressed through the investigation of the causal effects of the treatment histories 
A*{t- s + 1, t) on the outcomes Y* G L\t + 1) for dl\t&Tg. 

We want to underscore again that in this approach and for a given t E %, 
the outcome is not defined as a time-dependent variable in the sense that 
it corresponds with a variable measured at a unique time-point, specifically the 
last time-point t+1 associated with the corresponding t-spccific observed data, 
O*. That is why, although F* = Y{t + 1), we adopt a separate notation F* 
to designate the outcome. It is not to be confused with the notation Y intro- 
duced for the conventional counterfactual framework and which designates a 
time-dependent variable. Indeed, Y{j) can be regarded in the <-specific counter- 
factual framework both as any covariate Y{j) G L*{j) and the outcome F* when 
j = t + I. Similarly, note that we adopt a distinct notation, A* , to unambigu- 
ously represent the treatment of interest in the i-specific observed data 0{t). 
This notation is not to be confused with A which refer to a variable that can 
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be regarded both as a covariate i* and a treatment variable A* in the f-specific 

counterfactual framework. 

Note that if for example s = 2 , then 1^(1) cannot be of interest since at time 
point 1 each unit has been treated with a treatment history of size 1 and it is 
thus not possible to look at the effect of a treatment history of size 2 on Y{1). 
That is why we use the notation 7^ to indicate that the set of outcomes of in- 
terest depends on the investigator's choice for s. 



A. 2. Assumptions 

In the t-specific counterfactual framework, we adopt the same set of assumptions 
as described for the conventional counterfactual framework with the exception 
that each assumption is made relative to each t-specific observed data of inter- 
est, O* for t G 7^. In other words, we make the following t-specific assumptions 
for all tG%. 

Existence of counterfactuals: we assume the existence of the following t- 
specific treatment specific processes, £^t((_g_,_i ^^(t— s+1, t+1), also referred to as 
t-specific counterfactual processes, for every treatment regimen a*{t — s + l,t) ~ 
{a{t — s + 1), . . . ,a{t)) e Av{t — s + l,t) where Av{t — s + designates all 
possible treatment regimens between time points t—s+1 and t, i.e. the support of 
the conditional distribution of A{t—s+l, t) given A{t—s) and V, g{A{t—s+l, t) \ 
A{t — s), V^).We denote the so-called t-specific full data process associated with 
O* with X* = {V,L^,^^_^^-^^^^{t - s + l,t + l))at(t-s+i,t)eAv{t-s+i,t) and its 
distribution with Fx* ■ 

Note that the existence of the t-specific counterfactual process L^a^^^t-s+i i) (^^ 
s+1, t+1) for every treatment regimen a*(t — .s+1. t) G Av{t—s+l, t) implies the 
existence of the t-specific counterfactual processes -£gt((_g_|_]^ ^-j (t — s + 1, j '+ 1) = 

-^a'(t-s+ij),Ato+i),...,A'(t)(*~*+l' ^* every i = t-s+l, . . . , i-1 and 
every treatment regimen a' (t — .s + l, J ) = (a(t — s+1), . . . , a(7)) G Av{t — s + j) 
where Avit — s + designates all possible treatment regimens between 
time points t — s + 1 and j, i.e. the support of the conditional distribution 
of A{t-s + l,j) given A{t - s) and V, g{A{t - s + l,j) \ A{t - s), V). We have 
A.v{t-s+l,j) = {a{t-s+l,j) : 3a{t~s+l,t) G Av{t-s+l,t) a{t-s+l,j) = 
a {t — s + l,j)} for j = t — s + 1, . . . , t — 1 and Av{t — s + l,j) is thus entirely 
defined by Av{t—s+l, t). Similarly, Av{t—s+l, t) is entirely defined by Av{K). 

Consistency assumption: at any time point j such that t — s + l<j<t+l^ 

we assume the following link between the t-specific observed data and the t- 
spccific counterfactuals: i*(j) — i^t((_^^]^ (j)- Under this assumption, we 
have: O* = {A\t- s + l,t),L%i^t-s+i,t)i^- ^ + ^^* + ^)) = 'l>\A\t-s+l,t),X'), 
where ^* is a specified function of the t-specific full data process X*. This no- 
tation indicates that the problem can be treated as multiple (for each t G Tg) 
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missing data problems. Indeed, for each t (z only the i-specifie counterfactual 
associated with the observed treatment A*{t — s + l,t) is observed; the others 
are missing. 

Temporal Ordering assumption: at any time point j such that t — s + 1 < 
j < t + 1, we assume that any treatment-specific variable can only be affected by 
past treatments: iat(t_,+i^4)(j) = iat(t_s+i (i) for j = i - s + 1, . . . , t + 1, 
where L*a.*{t-s+i t-s)(^ — s + 1) = L*(i — s + 1). This assumption is typically 
implied by the data collection procedure: the covariate is measured after 

A\t- 1) and before A\t). 

Sequential Randomization Assumption (SRA): at any time point j such 
that t — s-|-l<j<t-|-l, we assume that the t-specific observed treatment is 
independent of the t-specific full data given the t-specific data observed up to 
time point j: A^j) ± X* \ A^t - s + l,j - l),L\t - s + l,j). Under the SRA, 
the t-specific treatment mechanism, i.e. the conditional density or probability 
of A*(t - s-l- l,t) given X\ g{A\t - s + l,t) | X*), becomes: 

giA\t-,t) I X*) = f[ g{A'{j) \ A'{t-,j ~ 

t 

n 9{A'U)\A'{t-,j~l),L'{t-,j)), 

where t~ = t — s -I- 1. The SRA implies coarsening at random Q) and thus 
the t-specific likelihood of the t-specific observed data factorizes into two parts: 
a so-called Fx* and part. The Fx* part of the likelihood only depends on 
the t-specific full data process distribution, and the part of the likelihood 
only depends on the t-specific treatment mechanism. As a consequence of this 
factorization of the t-specific likelihood under the SRA, we now denote the 
distribution of the t-specific observed data with Pp ^ gt and the likelihood of 
O* is: 

^xt part part 

£(0')''=^/(L'(t-s+l)) llT=t-s+2 fiL'U)\L'{t-s+l,j-l),A*{t-,+l,j-l)) g(A'{t-s+l,t)\X') . 



In addition, we denote the set of conditional densities or probabilities that de- 
fine the Fxt part of the likelihood, except for /(i*(t — s -I- 1)) with Qf t- 



R. Neugebauer et al. /Causal inference with HRMSMs 



144 



A. 3. Equivalence between MSM parameters in the t-specific 

counterf actual framework and HRMSM parameters in the 
conventional counterf actual framework 

Like in the conventional counterfactual framework, causal effects can be rep- 
resented based on parameters defined by MSMs in the t-specific counterfac- 
tual framework. Indeed MSM approach can be applied to all i-specific observed 
data, O*. We refer to an MSM associated with a given t-specific full data as 
a t-specific MSM. These t-specific MSMs are t-specific full data models, i.e. 
model of Fx* , which define parameters based on a feature of the distribution 
of the following counterfactual outcomes: Y^ti^i_i^j^i ly Typically and specifi- 
cally in this article, one is interested in average causal effects per stratum V* 
of the population which can be represented by causal parameters defined by 
MSMs for -Ej^^t (Yg*t(-j_g^]^ | V^) for t G 7^. We denote a causal parameter 
defined by such an MSM with l3t{Fxt \ ■) to indicate that it is a mapping 
from the space of t-specific full data distribution Fx* to the space of real num- 
bers and that this mapping is a function of modeling assumptions represented 
by. 

We have by definition from V^* = V{t — s + 1) and we can show as follows 
that = Fa(t-.+i,t)(i + 1) for t £ %: 

Y-aHt-s+i,t) = from © 

= Ya(t-s+i,t) from (P 

Thus we have Ep^, {Y-^^i^t-s+i,t) I = + 1) I V{t - s + 1)) 

and Fx* = ilj{Fx,g) for some specified function tjj. 

In general, one can show that HRMSM parameters defined in the conventional 
counterfactual framework, f3t{Fx,g \ ■) (see section [^3]) . corresponds to MSM 
parameters defined in the t-specific counterfactual framework, PtiFx* \ •): 

/3tiFx,9\-)^HFx^ \-) (6) 

In addition, note that MSM parameters defined in the t-specific counterfactual 
framework are typically different from MSM parameters defined in the conven- 
tional counterfactual framework. Indeed the former are defined as a function of 
both Fx and g while the latter are defined as a function of Fx only. 

A. 4- Link between the conventional and t-specific counterfactual 
frameworks 

Figure [5] illustrates based on an example of a longitudinal study with short 
follow-up the link between the longitudinal data representation in the conven- 
tional counterfactual framework and its representation in the time-specific coun- 
terfactual framework. Note that in the conventional counterfactual framework 



Example of a longitudinal study with K=3 and T=[0,1.2,3} 




For aU t in T„ we have O' = ( L' (t-1) , A' (t-1) , L' (t), A' (t), L' (t+1) ) = ( A'(t-s+l,t), L' (t-s+l,t+l) ) 
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the data are approached as a single entity, O, in the sense that the treatment is 
defined once and for all as a history A[K) and the outcome is time-dependent, 
Y{t) £ L{t). On the other hand, in the t-specific counterfactual framework 
the data are viewed as layers of separate entities, O*. For each O*, the treat- 
ment and outcome of interest are redefined along with the baseline covariates 
(highlighted in yellow/light gray on figure [5|). Note that for each O*, the out- 
come is no longer time-dependent but correspond with the last outcome col- 
lected, Y*- e L^{t + 1), (highlighted in orange/dark gray on figure [S]) and the 
treatment history size is fixed to a user-specified value s = 2. In the conven- 
tional counterfactual framework, the investigator examines the effect of A{K) 
on Y{t) for all i G T based on MSMs for the full data associated with O 
whereas in the i-specific counterfactual framework, the investigator can exam- 
ine the effect of A\t - s + l,t) on F* for t e % based on MSMs for the 
t-specific full data associated with O* . Figure [5] illustrates how the t-specific 
counterfactual framework can be viewed as a collection of conventional coun- 
terfactual sub-frameworks with distinct definition of the outcome, treatment 
and baseline covariates. These conventional counterfactual sub-frameworks dif- 
fer from the conventional counterfactual framework in the sense that the treat- 
ment history is of size s K + 1 and the outcome of interest is no longer 
time-dependent. 

We have previously underscored that the i-specific and conventional coun- 
terfactual approaches lead to a different representation of the causal effect of 
A on. Y with MSMs. We have also shown the equivalence between the pa- 
rameters defined with MSMs in the i-specific counterfactual framework and 
HRMSM parameters defined in the conventional counterfactual framework and 
argued that these parameters may often provide a representation of causal ef- 
fects which are more relevant for public health research. We can thus now eas- 
ily address the issue of HRMSM parameter estimation. HRMSM estimators 
can indeed be derived from the application of the conventional MSM estima- 
tion methodologies to all the aforementioned conventional counterfactual sub- 
frameworks. 



A. 5. HRMSM estimators 

Under the assumptions presented earlier in section [A. 21 and from equality ([6|), 
the HRMSM parameters, /3t{Fx,g \ ■), can be identified and consistently esti- 
mated with the t-specific observed data and three estimators of the f-specific 
MSM parameters: the Inverse Probability of Treatment Weighte d (flil. fioh . 
the G-computation ([ll[ll,|l,[ll,[2i,[2i,[8l[^ and Double Robust HEIMI^ esti- 
mators. The implementation procedures for these three estimators of HRMSM 
parameters correspond with the procedures developed for MSM-based causal 
inference except that they are applied not to the observed data, O, with treat- 
ment, A, and time-dependent outcomes of interest, Y(t) for t G T, but to all 
t-specific observed data, O*, of interest, i.e. for < G 7^, with treatment A* and 
outcome F*. The consistency and efficiency properties of these three estimators 
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along with their implementation procedures have been thoroughly studied in 
the literature (see references above). 

Appendix B: Sufficient assumptions for consistent HRMSM 
estimation 

We now formally establish that the seemingly larger set of i-specific assumptions 
(see section IA.2P required to investigate causal effects in the t-specific counter- 
factual framework is implied by the set of assumptions (see section [2T2|) required 
to investigate causal effects in the conventional counterfactual framework. The 
important practical consequence of both of these results is that successful in- 
vestigation of causal effects with HRMSM parameters can be achieved based on 
the same model assumptions leading to successful investigation of causal effects 
with MSMs. Thus, the choice of HRMSM-based causal analysis over MSM-based 
causal analysis is only a matter of practical considerations (statistical power and 
computing issues) and above all subject-matter considerations (the relevance of 
the causal effect representation to public health research). 

Theorem B.l. We adopt the notations introduced previously for the conven- 
tional and t-specific counterfactual frameworks and in particular t~ = t — s -\- 1. 
Based on these notations we have: 

i. the assumption of existence of counterfactuals defined in the conventional 
counterfactual framework implies the t-specific assumptions of existence of 
counterfactuals defined in the t-specific counterfactual framework: 

y a{K) e Av{K) L^^K){K +1) =^ \/ t e {s - 1, . . . , K} V a\t- ,t) G 

Avit~,t) L^.(,_^,)(t-,i-Hl) 

ii. the consistency assumption in the conventional counterfactual framework 
implies the t-specific consistency assumptions in the t-specific counterfac- 
tual framework: 

yte{0,...,K + i} L{t) ^ LA^K){t) =^ yte{s-i,...,K} Vje 

{<-,..., t + i} i*(.?) = i!i*(t-,t)(j) 

Hi. the temporal ordering assumption in the conventional counterfactual frame- 
work implies the t-specific temporal ordering assumptions in the t-specific 
counterfactual framework: 

Vte{0,...,X + i} ia(K)(i) = ia(t-i)(i) =^ Vte{s-i,...,K} 



yje{t-,...,t + i} 4'(t-,o W = 
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iv. the SRA in the conventional counterf actual framework implies the t-specific 
SRAs in the t-specific counterfactual framework: 

yte{0,...,K} A{t) ±X\A{t~l),L{t) =^ yte{s-l,...,K} 

y J e{t-,...,t} A\j) ±x'\A\t-,j-i), L\r,j) 

Proof. For i e {s — 1, . . . , K} and j G . . . , f + 1} we have: 
iv. 



X' = [v,Ll,^,-^^^{t-,t + l)^ 

= (v,Ll^^^^,)(t-,t+\)\ fromO 

V J a(t- ,t)eAv(t- ,t) 

= {v, A{t - s) , L-,<^t- ^t) (t-), L-,(t- ,t) {t- + l,t+ 1)^ 
from ^ and ^ 

= (v,A{t-s),L-,^t-^t){t + i)^ 
= (v,A{t-s),(L^^t-^f^{t+l)) ) 

\ \ / a{t- ,t)£Av{t- ,t)/ 

X' - {V,A{t-s),Xl)vjheTeXl^(v,L-,(t-^t){t + l)\ (7) 

\ J a(t- ,t)eAv(t- ,t) 



a*(t-,t)G^v(t-,t) 
a{t- .t)eAvit- ,t) 

I 

a(t- ,t)£Av(t- ,t) 

a{t- ,t)GAv{t^ ,t) 
I 

a{t- ,t)£Av{t- ,t)/ 



In addition, we have: 



= ^A{t-s),a{t- ,t),A{t+l,K}it + 1)^ 



a{t- ,t)£Av(t~ ,t) 



a{t- ,t)<£Av{t~ ,t) 



C X = [V,La(^K)iK+l) 



a{K)eAv{K) 

Xl c X (8) 
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Based on these previous two results, we obtain: 

gl^A\j),X'\A\t-,j~l)rL\t-^j)^ 

= gl^A{j),X' \ A{t- ,j - l),L\t-),L\t- + 1,])^ from O 
= g{A{j),X'\A{t-,j^l),A{t~s),L{t-),L{t- + from da) and 

= .g(^A(j),X*|A(j-l),Z(j)) 
- 9l^A{j),A{t - s),Xl I A{j - l),Z(j)) from 
= g(^A{j), Xl I A{j - 1), Z(j)^ since A{t - s) C A{j - 1) 
= g(^A{j) I A{j - from the SRA and ^ 

= g(A\j) I A\t-,j - l),L\t-,j)] from (P), © and © 



This last equality is equivalent to A*{j) ± X* \ j - l),L\t~,j). 



We also have: 
• iij^t-. 
i. 



4'(t-,t)(j) = 4(t-,t)(j) from o 
= ia(t-,t)0') from © 

= LA(0),...,A{t-s),a{t- A),A{t+l),...,A{K){j) 



L\j) = L{j) from © 

= La(k){3) fr'om the consistency assumption 

= -^A(t-s),^(t-,t),A(t+l,i<-)(j) 

= ^A{t-,t}{j) 

= LAt(t-,t){i) from © 
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4'(t-,t)0') = L-a(t-,t){i) from © and © 

= I^A(t-s),a(t- ,t),A(t + l,K){i) 

= iA(t-s),a(t-,i-i)0) from the temporal ordering 
assumption 

= 4'(t-,,-i)(j) from © and © 



^a'(t-,t)(j) = 4(t-,t)(j) fr-om © 



A{0,t- s),L{0,r)\ from© 

/ a{t~ ,t) 

(^A{t — s), iyi(o),...,yl(t-s),a(t-,t),A(t+l),...,yl(if)(i )^ 



L\j) = ( A(0,i-s),Z(0,<-)j from ® 

A{0, t — s), Z^(^-,(0, t^)^ from the consistency 
-ssumption 

A{0,t - s),iA(t-s),A(t-,i),^(*+i,Ar)('^'^ )^ 
(^A(0,t-s),Z^(,-,,)(0,i-)^ 
A(0,i-s),Z(0,<") 



A{t-,t) 



LAit-.t)ij) from® 
L%it-,t)(j) from © 
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4*(t",t)(i) = [A{0,t~ s),L{Q,t-)] from (HI and ® 

a{t- d) 



^(0,t-s),La(t-,t)(0,t-) 

^(0,t — s), Z^((_5-)(0, i^)^ from the temporal 
ordering assumption 
= {A{Q,t-s),L{Q,t-)\ 

V J A{t-s) 

= L*j(,„,)(j)from(l3l) 

= ^Ait~s),aHt-,j~i)ij) a\r,t- s) is empty by 

definition 

= 4*(t-j-i)(j)° 

Theorem B.2. FFe adopt the notations introduced previously for the conven- 
tional and t-specific counterfactual frameworks and in particular t^ = t — s + I. 
Based on these notations we have: 

I. \f te{s-_l,...,K} W j e{t-,...,t+i} 

f{L\j) I Z*(t-,j - l),A\t-,j - 1)) = f{L{j) I L{j - 1), A(j - 1)) 

a. \J t^{s~i,...,K} e {t- t} 

g{A\j) I A\t-,j - l),L\t-,j))^g{A{j) \ A{j - l),L(j)) 
Proof. 

For t E {s — 1, . . . , K} we liave 
i. for j e{t- + l,...,t+l}: 

f(^L\j)\L\t-,j^l),A\t-,j-l)^ 

= f(^L\j) \ L'{t-),L\t- + l,j ^ l),A{t- ,3 - 1)^ from© 

= f (^L{j) \ L\t-),L{t- + l,j - l),A{t- ,j ~ 1)^ from© 

= f(^L{j) I A{t-s),L{t-),L{t- + l),A{t-, J- 1)) from ® 

= /fit?) |ZO--l),A(j-l)' 
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for i : 

/(^L*(j) |Z*(t-,j-l),i*r,j-l)^ 

= f(L\j)\L'(t-),L\t- ^\,j^l),A{t-,j-l)^ from© 

- f{L\])\L'(t-),l{t~ + -I), A{t-, J -l)^ from© 

= fl^L\j)\A{t~s),L{t-),L{t- + l,3-l),A{t-,3-l)^ from© 

- / [a{j - ^),m I L{j - 1), ^(j - 1)) from m 

- f(Lij)\L{j~l),Aij-l) 



ii. for j e {t , . . . ,t}: 

5(^*0-) i^*r,j-i),^*r,j)) 

= M(t-,j-l),L*r),Z*(t- + l,j)^ from© 

= 5(^0) + from© 

= g(^A{j)\ A{t-, J- 1), A{t-s),L{r),L{t- + from® 

= g[A{j)\A{j-l).L{j) 
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